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Abstract 

We describe the software package SPEX, which allows first-principles cal- 
culations of quasiparticle and collective electronic excitations in solids us- 
ing techniques from many-body perturbation theory. The implementation 
is based on the full-potential linearized augmented-plane-wave (FLAPW) 
method, which treats core and valence electrons on an equal footing and 
can be applied to a wide range of materials, including transition metals 
and rare earths. After a discussion of essential features that contribute 
to the high numerical efficiency of the code, we present illustrative results 
for quasiparticle band structures calculated within the GW approximation 
for the electronic self-energy, electron-energy-loss spectra with inter- and 
intraband transitions as well as local-field effects, and spin-wave spectra 
of itinerant ferromagnets. In all cases the inclusion of many-body correla- 
tion terms leads to very good quantitative agreement with experimental 
spectroscopies. 

1 Introduction 

First-principles computations, which do not rely on any empirical input pa- 
rameters, have become an important tool in materials science. Ideally, such 
calculations use only the fundamental laws of nature together with the specified 
elemental composition in order to predict the structural and electronic proper- 
ties of a material. As chemical bonding and the response to external fields are 
determined by the microscopic dynamics of the ions and electrons inside the 
solid, the relevant laws in this case are those of quantum mechanics. The com- 
plete Hamiltonian can, in fact, be readily written down, because the interaction 
potentials between all constituent particles, given by Coulomb's law, are known 
exactly. However, a direct solution of the Schrodinger equation or its relativistic 
counterpart, the Dirac equation, is not feasible for extended solids, because the 
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number of electrons iV is of the order of and hence too large for a nu- 

merical treatment of the correlated many-electron wave function Vl/(ri, . . . , rjv). 
Instead, practical applications rely on alternative approaches that are formally 
equivalent to the Schrodinger equation but do not employ the many-electron 
wave function as the basic variable. A prominent example is density-functional 
theory (DFT), which is based on the ground-state electron density n(r), a real 
quantity that depends only on a single spatial position vector. In spite of this 
vast simplification, the Hohenberg-Kohn theorem [T] asserts that knowledge of 
the ground-state density alone is sufficient, at least in principle, to determine 
all observables of a stationary system in equilibrium. The density itself can be 
calculated from the Kohn-Sham scheme [2], which introduces an auxiliary sys- 
tem of non-interacting electrons with the same ground-state density as the real 
interacting system and requires the solution of a single-particle Schrodinger-like 
equation, the Kohn-Sham equation, with a self-consistent local potential. 

Despite the enormous scope of the Hohenberg-Kohn theorem, practical ap- 
plications are limited because for most observables no explicit formulas for the 
actual dependence on n(r) are known. A notable exception, besides the den- 
sity itself, is the ground-state total energy, where the dominant contributions 
are given exactly and the remaining exchange-correlation functional may be re- 
placed by the local-density approximation (LDA) |2j or generalized gradient ap- 
proximations like PBE [3] . By comparing the total energies for different atomic 
configurations it is thus possible to predict crystal structures and related quan- 
tities like the elastic moduli without empirical parameters. The huge success of 
density-functional theory stems from the good agreement with crystallographic 
measurements for a wide range of different materials. In contrast, electronic 
excitation spectra often show substantial deviations from experimental spectro- 
scopies. The most famous example is the severe underestimation of the band 
gaps of semiconductors. As a large part of the error in this case is systematic 
and arises from the pervasive but incorrect identification of the Kohn-Sham 
eigenvalues with the true quasiparticle energies, better approximations for the 
exchange-correlation functional will not solve the problem. Therefore, quanti- 
tative methods for electronic excitations in solids are now chiefly sought outside 
density-functional theory. Perhaps the simplest extension are hybrid functionals 
like PBEh [4] or HSE [5] , which replace a fraction of the local exchange potential 
by non-local Hartree-Fock exchange. As Hartree-Fock in turn tends to overes- 
timate band gaps, a linear combination with suitably chosen coefficients can 
improve the agreement with experiments. Unfortunately, the optimal relative 
weights depend on the material: While an admixture of 25% Hartree-Fock yields 
good results for many typical semiconductors, a higher fraction is required for 
large-gap insulators, where screening is much weaker [5]. For metals any inclu- 
sion of Hartree-Fock exchange even increases the already too large valence-band 
widths further, so that hybrid functionals are in fact in worse agreement with 
experimental band structures than pure LDA results |6J. 

A better-founded scheme without adjustable parameters is many-body per- 
turbation theory [7], which is based on the Green function G(r, r';cj). As it 
contains no inbuilt systematic errors, the accuracy is only limited by functional 
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approximations like the GW approximation for the electronic self-energy [S], 
which can be improved if necessary. As a consequence, this approach yields exci- 
tation energies in significantly better agreement with experiments than density- 
functional theory. The main drawback is the high computational cost, which 
is related to the more complicated mathematical form of the non-local and 
frequency-dependent Green function in comparison to the density. For this 
reason, practical applications have so far been limited to moderately complex 
systems with up to about one hundred atoms per unit cell. Although this suffices 
to study, for example, carbon nanotubes [5] or point defects at semiconductor 
surfaces flOj. simulations with tens of thousands of atoms are now possible us- 
ing density- functional theory with linear-scaling algorithms Furthermore, 
many actual GW calculations contain a number of additional simplifications. 

Besides individual quasiparticle properties, collective excitations of the elec- 
tron system in solids constitute another major challenge for electronic-structure 
calculations. Prominent examples are plasmons and excitons, which are associ- 
ated with resonances in the dielectric function and can be probed by electron- 
energy-loss or optical spectroscopies, but also spin-wave modes in magnetic ma- 
terials. As a single-particle picture cannot describe such collective excitations, 
many-body perturbation theory has become the method of choice for quantita- 
tive simulations \12\ . 

Here we describe a new software package, SPEX [T3], which contains an im- 
plementation of many-body perturbation theory and can be used to simulate 
various spectroscopic techniques. The code is designed to avoid many short- 
comings of previous implementations and to keep additional approximations at 
a minimum. Most importantly, it is not based on the prevalent pseudopotential 
concept but uses the full-potential linearized augmented-plane-wave (FLAPW) 
method [TJ, which opens the door to a wider range of materials, including 
elements with localized d or / orbitals. In the following we first discuss some es- 
sential features of the code. Then we present selected results for several spectro- 
scopies that probe electronic excitations, such as quasiparticle band structures, 
electron energy loss or spin waves, in order to illustrate the possible applications. 
Finally, we briefly summarize and comment on planned future developments. 

2 The SPEX code 

Although many-body perturbation theory is a self-contained mathematical frame- 
work, actual applications are most efficient when combined with density-func- 
tional theory. The Kohn-Sham eigenvalues, which are often in qualitative agree- 
ment with the true band structure, then serve as the starting point for a dia- 
grammatic expansion in terms of the dynamically screened Coulomb interaction 
W{r,r';uj). The screening arises from the formation of exchange-correlation 
holes around charged fermions and leads to the concept of quasiparticles, which 
incorporate the polarization of the local environment and constitute one class 
of elementary excitations of the electron system. As the screened interaction 
between the quasiparticles is much weaker than the bare Coulomb potential, a 
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Figure 1: Flowchart for the perturbative calculation of different electronic exci- 
tations and spectroscopies as implemented in SPEX. 

perturbative treatment is justified. A first-order expansion of the self-energy 
leads to the GW approximation, which is just the linear term in W{r,r';Lu). 
In accordance with standard perturbation theory, its matrix elements must be 
evaluated with the original unperturbed Kohn-Sham orbitals. 

Following this philosophy, SPEX is designed as a separate module that com- 
putes electronic excitation spectra from a given set of Kohn-Sham orbitals and 
eigenvalues. It can be combined with any density-functional code whose output 
data is convertible to the FLAPW form. We use our own package FLEUR jl3j . 
but other choices are equally possible. The flowchart in Fig. [T] illustrates the 
course of the calculation. 

A premier goal during the code development was to avoid additional sim- 
plifications wherever possible, so that the results depend exclusively on con- 
trollable convergence parameters and the choice of functional approximations 
like the GW approximation for the self-energy or the random-phase approx- 
imation (RPA) for the dielectric function. Besides, emphasis was placed on 
high efficiency, especially in terms of CPU time, so that complex materials with 
large unit cells can be treated. Finally, the code should be versatile and eas- 



4 



ily adaptable to new physical problems, for instance in the emerging areas of 
nanomagnetism and spintronics. To satisfy these requirements SPEX contains 
a number of important features that are summarized below. 

(i) Consistent employment of FLAPW, which treats core and valence elec- 
trons on the same footing. Unlike the atomic-sphere approximation used in 
early all-electron codes, the full-potential treatment is also suitable for surfaces 
or defects. The linearization error caused by expanding the muffin-tin func- 
tions around fixed energy parameters, which is negligible for states close to the 
Fermi level but becomes important for high unoccupied bands that contribute 
to the self-energy, can be eliminated by including higher energy derivatives as 
additional local orbitals |15| . 

(ii) The FLAPW basis set is optimized for the Kohn-Sham orbitals but does 
not span the complete Hilbert space. Representing products of wave functions, 
for example in the polarizability or matrix elements of the Coulomb potential, in 
this basis hence implies a loss of accuracy. Instead, we use a mixed product basis 
[TB] that is constructed from products of the basis functions and allows exact 
projections. After linear dependencies are removed, this set may be truncated 
in a controlled fashion, if desired. 

(iii) The dielectric function is constructed either in the random-phase or the 
time-dependent local-density approximation, without recourse to plasmon-pole 
models. The full frequency dependence as well as local-field effects are thus 
properly included. Besides, it gives access to the imaginary part of the self- 
energy, which leads to complex quasiparticle energies and describes the finite 
lifetime of the excited states due to scattering. 

(iv) Frequency convolutions are normally evaluated by means of contour in- 
tegrations in the complex plane, although the analytic continuation of functions 
calculated on the imaginary axis to real frequencies, which is faster but less well 
controlled, is also possible. 

(v) Symmorphic and non-symmorphic spatial symmetries are exploited wher- 
ever possible. For systems with inversion symmetry we also use the fact that 
the Coulomb matrix and response functions on the imaginary frequency axis 
may be chosen real and thus processed in compact form. 

(vi) Following a procedure developed in Ref. [17J for plane waves, the sin- 
gularity of the Coulomb matrix at k = in reciprocal space is treated exactly 
through an expansion that separates the divergent and regular parts [18| . A sub- 
sequent diagonalization confines the singularity to a single divergent eigenvalue, 
which can be processed analytically. At the same time, the transformation from 
the mixed product basis to the eigenvectors of the Coulomb matrix allows a 
very efficient truncation by eliminating the least significant scattering channels 
with the smallest eigenvalues. 

(vii) The spin degree of freedom is fully supported. We make no simplifying 
assumptions about the orbitals in the two spin channels and allow a completely 
spin-unrestricted treatment of magnetic materials. 

(viii) The code is applicable both to insulators and to metals. Where intra- 
band transitions require a special treatment, appropriate provisions are made. 
In particular, a Drude term is included in the dielectric function in this case. 
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All calculations are done at zero temperature. 

(ix) Relativistic corrections are treated at the scalar-relativistic level for 
valence states and the full Dirac equation for core states. 

As a performance test we conducted GW band-structure calculations for di- 
amond supercells of various sizes, using a k-point sampling equivalent to 4x4x4 
mesh points in the full Brillouin zone corresponding to the elementary diatomic 
unit cell and with optimized parameters that otherwise guarantee a convergence 
of the quasiparticle shifts to within 0.01 eV 112]. Our results demonstrate that 
simulations even with 128 atoms per cell are perfectly feasible on a standard 
single-processor work station, and a complete quasiparticle band structure re- 
quires less than one and a half days of CPU time in this case. We find a scaling 
behavior that lies between quadratic and cubic with the number of atoms in 
this size range. 

3 Results 

In this section we present illustrative results for different spectroscopies. All 
materials discussed here contain transition-metal elements and would be difh- 
cult to treat with conventional pseudopotential codes, because the localized d 
orbitals require a large number of plane waves. 

3.1 Quasiparticle band structures 

In the perturbative approach adopted here the quasiparticle energies are derived 
from the Kohn-Sham eigenvalues and the self-energy according to 

6„k. = e^^L + {^'^L\^A^n^..) - KflV'J^kJ • (1) 

The self-energy is evaluated in the GW approximation, and the matrix elements 
of the local exchange-correlation potential V^'^(r) are subtracted to avoid double 
counting. Although GW calculations for real materials have been feasible since 
the 1980s |i20ji21j, the prevailing reliance on plane waves and pseudopotentials 
meant that applications were long restricted almost exclusively to sp-bonded 
semiconductors and simple metals. No such restrictions apply to the FLAPW 
method. As an example, in Fig. [2] we show the band structure of strontium 
titanate (SrTiOs), a high-K dielectric insulator from the perovskite family, in 
the cubic phase observed at room temperature. The lattice constant within 
DFT-PBE [3j is 7.46 Bohr. For the self-energy construction we use 64 mesh 
points in the full Brillouin zone, 550 unoccupied bands, and cutoff parameters 
Gjnax = 5 Bohr" ^ and £max = 4 for the mixed product basis in the interstitial 
region and in the muffin-tin spheres, respectively. Second energy derivatives 
are included as additional local orbitals. Although DFT-PBE gives the correct 
qualitative picture, both the indirect band gap of 1.81 eV between R and F 
and the direct gap of 2.18 eV at F are significantly too low. In contrast, the 
GW approximation yields 3.23 eV and 3.61 eV in very good agreement with 
the experimental values 3.25 eV and 3.75 eV [22]. The size of the band gap 
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Figure 2: Quasiparticle band structure of cubic SrTiOa calculated in DFT-PBE 
(dashed lines) and the GW approximation (solid lines). 



is quite important in this SrTiOa has long been technically used as 

an optically transparent synthetic diamond simulant. The comparison with a 
previously published GW value of 5.07 eV for the indirect gap [53], obtained with 
a parametrized model dielectric function, underlines the need for an accurate 
evaluation of the self-energy correction. 

In Fig. [3] we display results for a wider range of materials from small-gap 
to large-gap insulators. In all cases the GW approximation corrects virtually 
the entire error of the LDA. The data support our previous conclusion, based 
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Figure 3: Calculated band gaps of selected materials in the LDA (circles) and 
the GW approximation (squares) compared to experimental values. 
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on an in-depth study of silicon [TS], that carefully performed all-electron GW 
calculations yield very good agreement with experiments. This had initially 
been in doubt after an early FLAPW calculation found unexpectedly large de- 
viations both from experiments and from established pseudopotential values 
|24| . However, this observed discrepancy is now understood to have arisen from 
incomplete convergence, especially with respect to the number of unoccupied 
conduction bands in the self-energy (TSl [5S]. Our own results are meanwhile 
confirmed by other independent all-electron implementations |26l I27j . 

Incidentally, these studies also revealed that well converged all-electron GW 
band gaps do not coincide with pseudopotential values after all. The remaining 
difference stems from the pseudization of the wave functions and the linearized 
core- valence interaction in the latter approach [28) . 

3.2 Electron-energy- loss spectroscopy 

While the single-particle Green function and self-energy describe the final states 
of photoemission experiments, where the particle number changes due to elec- 
tron emission or injection, low-energy spectroscopies involving intraband or in- 
terband transitions are related to the dielectric function e(r, r',cj), which is 
linked to the density-density correlation function and characterizes the linear 
response to an external electric field. In particular, electron-energy-loss spec- 
troscopy (EELS) measures the imaginary part of the inverse macroscopic di- 
electric function, the so-called energy-loss function — Ime~^(q, a;). In Fig.|3]we 
display the EELS spectrum of ferromagnetic nickel calculated in the RPA. The 
response function is constructed with 40x40x40 mesh points in the full Bril- 




Figure 4: Electron-energy-loss spectrum (EELS) of ferromagnetic Ni for q = 
(0.25, 0, 0)27r/a calculated in the RPA with (solid line) and without (dot-dashed 
line) transitions from the 3s and 3p core states compared to experimental data 
(dashed line) from Ref. [29j. 
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louin zone, Gmax = 5Bohr~^, Lmax — 4 and 118 unoccupied bands [TS]. Second 
and third energy derivatives are included as local orbitals. The resulting curve 
is in very good agreement with experimental measurements [i29], especially if 
transitions from the 3s and 3p core orbitals are taken into account; in this case 
the step in the energy- loss function around 64 eV, which corresponds to the on- 
set of transitions from these states, is also well reproduced. The RPA formally 
corresponds to a complete neglect of dynamic exchange-correlation effects in 
the linear density response function. The latter can be approximately included 
by the adiabatic local-density approximation [3Dj , but in most cases the results 
change very little compared to the RPA. 

3.3 Spin-wave spectra 

Electric fields couple to the charge of the electrons and induce characteristic 
excitations, such as interband optical transitions or collective plasmon modes, 
which correspond to resonances in the dielectric function and can be measured 
with frequency-resolved spectroscopies. Likewise, magnetic fields couple to the 
spin of the electrons and give rise to excitations of the spin system. These 
also fall into two groups, spin-flip Stoner excitations of individual electrons 
and collective spin waves, and can be identified as resonances in the dynamic 
transverse spin susceptibility. 

Most theoretical studies of spin waves are based on the Heisenberg model, 
although the assumption of localized spins makes its justification doubtful for 
ferromagnetic metals with itinerant electrons. In contrast, very few attempts 
at first-principles calculations were reported so far [31', "35]. The difficulty 
is twofold: As magnetic behavior originates in localized d or / orbitals of 
transition-metal and rare-earth elements, an all-electron scheme is mandatory. 
Furthermore, the treatment of spin waves requires dynamic exchange-correlation 
effects not contained in the RPA. We have solved the latter by explicitly in- 
cluding the screened Coulomb interaction between electrons and holes in the 
two different spin channels |33) . The resulting two-particle problem is analo- 
gous to the Bethe-Salpeter equation used for excitons in semiconductors. As 
the dominant part of the dynamic correlation in ferromagnets is caused by mul- 
tiple scattering of electron-hole pairs at the same atomic site, we transform the 
screened interaction to a basis of maximally localized Wannier functions |34) . 
which allows a very efficient truncation. Here we take only matrix elements with 
four Wannier functions that are all localized at the same site into account, but 
systematic extensions are of course possible to ensure full convergence [35J. 

As an example we show results for iron, obtained with 30x30x30 mesh 
points in the full Brillouin zone, Gmax = 4.5Bohr^^, Lmax = 4 and 100 un- 
occupied bands. Figure [S] shows the transverse spin susceptibility of the non- 
interacting Kohn-Sham system. As there is no dynamic correlation in this case, 
only single-particle Stoner excitations exist. As a consequence, the spectral 
function ImX]^<t(q, oj) exhibits a peak at around 2eV, which equals the ex- 
change splitting visible in the density of states and corresponds to spin-flip 
transitions between occupied majority and unoccupied minority states. When 
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Figure 5: (a) Density of states DOS(cj) for the majority (up) and minor- 
ity (down) spin channel of ferromagnetic Fe. (b) Imaginary part of the non- 
renormahzed Kohn-Sham spin susceptibility for q = (0.1, 0.1, 0)27r/a. 



dynamic correlation is included in the form of the screened Coulomb interaction, 
an additional spin-wave peak appears at low energies as illustrated in Fig.iJa). 
Plotting the peak positions as a function of the wave vector yields the spin- wave 
energies as displayed in Fig. IHl^b) for the [110] direction in iron. The dispersion 
obtained in this way, with no empirical parameters, is in excellent agreement 
with experimental data ^36j. 




CO (eV) 

Figure 6: (a) Imaginary part of the renormalized spin susceptibility of Fe for 
wave vectors q = (^, ^, 0)27r/a and (b) calculated spin- wave dispersion along the 
[110] direction compared to experimental data from Ref. |36) . 
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4 Summary and outlook 



We have discussed a new implementation of many-body perturbation theory 
within the FLAPW method. A number of features, such as the copious use of 
symmetries, the treatment of the Coulomb singularity or expedient basis trans- 
formations that allow efficient truncations, which are described in detail else- 
where jl8! , ensure a high computational efhciency. Our results for electronic 
excitations and associated spectroscopies in solids are in excellent quantitative 
agreement with experiments and show that full-potential calculations are now 
feasible even without plasmon-pole models or other far-reaching simplifications. 
As the FLAPW method is applicable to transition metals and other complex ma- 
terials that cannot be easily treated with pseudopotentials, this opens up new 
prospects for theoretical investigations. We have already started out on this 
path by exploring spin-wave excitations in ferromagnets. Without any empiri- 
cal parameters, we obtained spin-wave dispersions in very good agreement with 
experiments. The rich physics of spin-dependent phenomena means that many 
further developments are still necessary, however. For example, the relativistic 
spin-orbit coupling and non-collinear magnetism require appropriate extensions, 
as do spin dynamics at finite temperature or the inclusion of spin-dependent 
scattering in the electronic self-energy. Outside the field of nanomagnetism, 
the linear and non-linear optical properties of modern optoelectronic materials 
provide another timely challenge. 
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